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We study the susceptibility propagation, a message-passing algorithm to compute 
correlation functions. It is applied to constraint satisfaction problems and its accu- 
racy is examined. As a heuristic method to find a satisfying assignment, we propose 
susceptibility-guided decimation where correlations among the variables play an im- 
portant role. We apply this novel decimation to locked occupation problems, a class 
of hard constraint satisfaction problems exhibited recently. It is shown that the 
present method performs better than the standard belief-guided decimation. 

I. INTRODUCTION 

Message-passing algorithms have shown to be effective in helping to find solutions of 
some hard constraint satisfaction problems (CSPs) like if -satisfiability and coloring. The 
simplest application consists in using belief propagation (BP), when it converges, in order 
to get some estimate of the marginals of each of the variables. One must then exploit the 
information obtained in this way (which is in general only approximate in a CSP described by 
a loopy factor graph). So far only two methods have been explored thoroughly: decimation 
and reinforcement [3j]. Decimation consists in identifying from some criterion the most 
"polarized" variable (e.g. the one with the smallest entropy), and in fixing it to its most 
probable value. After this variable has been fixed, one obtains a new, smaller, CSP, to 
which one can apply recursively the whole procedure (BP followed by identifying and fixing 
the most polarized variable). In reinforcement, one finds from the BP marginals the most 
probable value of each variable, and one adds, in the local measure of each variable, an extra 
bias in this preferred direction. The new CSP therefore has the same number of variables 
as the original one, but the local measure on each variable has been changed. One iterates 
this reinforcement procedure until the variables are infinitely polarized. If the algorithm 



is successful this returns a configuration of variables which satisfies all constraints. These 
two procedures, BP+decimation and BP+reinforcement, are remarkably efficient in random 
CSPs like /^-satisfiability [2], graph colouring |2|, and perceptron learning |4]. When one 
approaches the SAT-UNSAT threshold of these problems, a more elaborate version which 

nnrt 

uses the information on marginals from survey propagation (SP) is more effectivepj, |3|, |5j, 
and at present the SP-based decimation and reinforcement methods are the most efficient 
incomplete SAT solvers for random 3-satisfiability. 

Recently, a class of problems has been described jf| [?| where these procedures are much 
less efficient. These are the locked occupation problems(LOPs), a class of CSPs where the 
set of solution consists of isolated configurations, far away from each other. Apart from the 
XORSAT problem^ which can be solved by Gaussian elimination, the random LOPs are 
very hard to solve in a broad region of the density of constraints, below their SAT-UNSAT 
transition. For these LOPs, it is known that SP is equivalent with BP. The BP+decimation 
method has been found to give rather poor results, and the BP+reinforcement, which works 
better, is still rather limited. One reason for this hardness is the fact that local marginals 
often convey little information on the solution. This has motivated us to explore some 
extensions of the message-passing approaches, in which one uses, on top of local marginals, 
some correlation properties of the variables. Several possibilities to obtain information on 
the correlations from message-passing procedures have been explored recently 
Here we use the susceptibility propagation initially introduced in J4]. We show that some 
of the hard LOPs that could not be solved by previous methods can now be solved by a 
mixture of the single- variable decimation with a new pair-decimation procedure which makes 
use of the knowledge of correlation. In the case of binary variables which we study here, 
this new procedure amounts to identifying a strongly correlated pair of variables, and fixing 
the relative orientation of the two variables. 

The paper is organised as follows. In Section [Til we introduce the susceptibility propaga- 
tion, derived as a linear response to belief propagation. This method is examined analytically 
in Section IITT1 where it is applied to simple systems for which exact fixed points of the itera- 
tion are determined. In Section IIVI it is applied numerically to locked occupation problems 
and the accuracy of the method is examined: we measures the performance of the decimation 
process which makes use of the correlations obtained with this method. The final Section [V] 
is devoted to conclusion and discussions. 
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II. SUSCEPTIBILITY PROPAGATION 



A. Occupation Problems 

Let us consider an occupation problem, which consists of \V\ = N binary variables 
Xi G {0, 1} (i G V) and|F| = M constraints tp a (x ial , . . . , Xi a k ) — 1 (a G F). Each constraint 
involves exactly k variables and is parameterized by a (& + l)-component "constraint- vector" 
A = {A(0), . . . , A(k)) with binary entries defined as follows. We say a variable is occupied 
if Xi = 1. Let r a = ^2 ieda Xi be the number of occupied variables that are involved in the 
constraint ip a . By definition, the constraint a is satisfied (ip a = 1) if and only if A(r n ) = 1. 

An occupation problem is locked if the following three conditions are met ((]] {?)] jl2j] 

• A{0) = A{k) = 0. 

• A(r)A(r + 1) = for r = 0, . . . , k - 1. 

• Each variable appears in at least two constraints. 

Standard examples of locked occupation problems include positive 1-in-K satisfiability 
[13I 1 and parity checks [jjj]. 

As can be done for general constraint satisfaction problem, a factor graph G = (V, F; E) 
can be associated with an instance of the occupation problemsjisl]. The set of vertices 
of this bipartite graph G is V and F while the set of edges is E = {(i,a)\i G V, a G 
F, Xi is involved in ip a }- The notion of neighborhood is naturally introduced: da = {i G 
F\(i,a) G E}, di = {a G V\(i,a) G E}. For a collection of variables in S C V, we shall 
write x s = {xi\i G S}. We also use the short-hand notation x = x v . 

B. Belief Propagation Update Rules 

Consider an occupation problem described by a factor graph G = (V, F, E) and a 
constraint- vector A. For later use, we introduce local 'external fields' ft| (x G {0, 1}, i G V), 
which will be sent to zero at the end, and consider a joint probability distribution 

M N 
^— ' a=l 1=1 x 



This probability distribution is well defined as soon as there exists at least one ("SAT") 
configuration satisfying all the constraints. The constant Z(h x ) is a normalization factor. 
Our final aim is to extract solutions from the uniform measure p(x\0) over solutions satisfying 
all constraints (when there exists at least one solution). 

The marginal distribution pi(xi\h) can be estimated by the BP algorithm. The BP update 
rules for two families of messages, namely cavity fields and cavity biases, are given by 16.[l7] 

v£?(*i\k*) =^Ar^ n ^iMm x n ( 2 ) 

^i-^a\UL ) b€di\a x 

^m=J2 6 ^'M^ a ) n 4l a (^m- (3) 

x'g a idda\i 

Here, we have decided to introduce a normalization factor Z^ a (h x ) for v\% a (xi\lf) an d to 
avoid the normalization for v^^ix^h?) . This choice is perfectly valid for BP, and it helps 
to get relatively simple susceptibility propagation update rules ©©• 

Assuming convergence to a fixed point, the BP estimate for the marginal distribution of 
variable % is: 



n — ' bedi 



where v^i( x i\h?) ls t ne fixed point of the BP iteration. 



C. Susceptibility Propagation Update Rules 



The 2-point connected correlation function at h = is obtained as 

dpi(xi\h x 



Pij — Pij(Xi,Xj) Pi^Xj^jPjiXj) 



dhj 



(5) 



h=0 



To have a message-passing algorithm to calculate this quantity, we introduce the cavity 
susceptibility and its companion by 



Va~+i,j {Xi i Xj ) 



dvi^ a (xi\ff 



dv a ^i(xi\h x ) 



dhj 



h=0 



h=0 



(6) 



(7) 



Note that the roles of variables Xj and Xj are asymmetric; j can be an arbitrary variable 
while i is a neighbor of the constraint a. 



The cavity susceptibility and its companion can be calculated by a message-passing 



method 



9|. The susceptibility propagation update rules can be obtained by differentiat- 



ing the belief propagation update rules (j2J) and (j3J) with respect to h*. They read 4j||l8l] 



(t+i) / 



7 W (h x \ J - J - I ij (t 

^i^aVLL ) b£di\a \ b£di\a "b^iK^i) 

t&da\i I m(zda\i 



x' ga \e&da\i J meda\i v m^a\X 

where 

vt<k*i) = 4L{*i\k x = o), P^fa) = 4lM\h x = o). (io) 

The function C^_l a j(xj) originates from the derivative of z\% a and can be determined by 
requiring the normalization 

E^>^) = °- (n) 

Let us suppose that we have found a fixed point of BP and the susceptibility propagation. 
By differentiating (j3J) with respect to the external fields, we can express the 2-point connected 
correlation function in terms of the messages at the fixed point as 

b£di cedi\b 

The constant Cij(xj) is related to the derivative of Zi(h) and is conveniently fixed by the 
condition J2 X .pfj nn (xi, Xj) = 0. 



D. Log-likelihood representation 

The rules ( IH1I91) apply to all types of CSPs with discrete variables. When dealing with 
binary variables, it is helpful to rewrite the belief and susceptibility update equations in 
terms of log-likelihood variables. We introduce the cavity field and cavity bias in the log- 
likelihood representation rii^ a and n a ^i as (we omit the time superscript (t) where it is 
obvious) : 

u i ^ a (x i \h)=A-, a e ni ^ s % (13) 
t>a^i(xi\h) =B a ^ e^*®* (14) 



where Sj is the spin variable Sj = 2xj — 1 = ±1 and the external fields in the two represen- 
tations are related by 

h 1 - h° 
j ~ 2 ' 

Naturally we define the cavity susceptibility in the log-likelihood representation as 



dni^ a (h) 



) 'la— n,j 



h=0 



dh,-, 



h=0 



The belief propagation update rules read 



n 



(t+i) 



=/a-*i({n,-_ +0 }j e 9o\i), 



(15) 



(16) 
(17) 



where 



fa-H({nj->a}jeda\i) = n 1° 



K-l) 



j£da\i 



By differentiating both sides of (I16fl7p . we obtain 



(t+i) 



pit) 



j£da\i) (t) 



m£da\i 



dn r 



(18) 
(19) 



(20) 
(21) 



Assuming that a solution n^ a of the BP equations ( TTBIfTTj) is used, one sees that the sus- 
ceptibility propagation update rule (I20|l20p is an inhomogeneous linear system in r\ and fj. 
The coefficient matrix takes the following form: 



dn r 



(22) 



1 " <^> 2 

where i,m G da and (•) means Here (sj) and (s m Si) for i,m G da means the expectation 
value with respect to the joint probability distribution for variables that are neighbors of a 
constraint obtained solely from beliefs[17|, Sec. 14.2.3]. 



In the log-likelihood representation, the magnetization and the pair correlation are given 
in terms of the fixed-point messages by 



fa) =tanh I J^n^ ) , 
1-tanh 2 



(23) 
(24) 



In the above expression, i and j can be arbitrary variables on the factor graph. 

III. PROPERTIES 

A. Linear Equation 

In order to study the structure of susceptibility propagation update rules ( l20f2Tl) . we 
construct a fcMiV-component column vector 



M _ (Jf) Jt) \t 



y'li^aj > 'la—>i,j)(i,a)eE,jeV • 



(25) 



Then the fixed point condition associated with ( l20f2~Ti) can be written as a linear equation 



y W = My 1 ** + b, 



with the inhomogeneous term 



b — (5ij, 0)(i t a)eE,jev 
The coefficient matrix is block- diagonal in j: 



M 



(iaj),(i'a'j') — 5j,j'-M-(ia),(i'a')j 



( 



M 







l(a' G di \ a)8i^ 




(26) 
(27) 

(28) 
(29) 



where the block A4 is independent of the block index j. 
Thus we obtain the unique fixed point 



y W = (1 -M^b 



(30) 



if (1 — M) is invertible, which is equivalent to the invertibility of (1 — M). 



The susceptibility propagation update rules (I20f21l) can be regarded as an iterative 
method to solve the linear equation equation flHUj) . It converges to a value irrespective 
of the initial vector if all the eigenvalues of M. have moduli smaller than unity. Because the 
block M does not depend on j, the existence of the fixed points and convergence to them 
are solely determined by M. and do not depend on j. 



B. Application to simple problems 

When the factor graph is a tree, even in presence of the external fields h x , the exact 
marginals are obtained by on a fixed point z^ a , z^2 a [l5(. Therefore, by differentiation 
with respect to h x , there exists a susceptibility fixed-point which gives the exact 2-point 
correlation function. In the examples which we have considered, the iteration of susceptibility 
propagation converges to this fixed-point. On the other hand, if the graph has more than 
one loop, there is no guarantee either that the fixed point exists or the iteration leads to 
that fixed point. In order to test these statements, we have studied a simple problem, the 
l-in-2 satisfiability problem, or anti-ferromagnetic Ising model. 

We first study this problem on a chain of length N. Namely, we take k = 2 and A = 
(0, 1, 0), and V = {0, 1, 2, . . . , N- 1}, F = {0 + ±, 1 + §, . . . , N- 1 - §}. This gives a simple 
case of a tree factor graph with E = {(i, = 0, 1, . . . , N— 2}U{(z, = 1, . . . , N—\}. 

Away from the boundaries,, since da and di consist of only two variables and constraints, 
respectively, (IToD, (fTTjl , (|2"DI). pT]) are simplified to yield 



i it+1 \ = n {t \ ; h {t \ =-n (i) 

i^i±7; t^^-ii i±7^—*i i±l— >i± 



ftW =-n«> (31) 

n {t+1 \. = t\ ,+Sij ; t\ .. = -<> (32) 

On the boundary, on the other hand, one has 

n 1 = n 3 = 0, 77 1 = Son, Ti 3= 5,- AT-l- (33) 

This in turn implies that 

J " (34) 
(otherwise) 



which gives: 



(si) = , (siSj) 



j /conn 



(*) (*) 

-n i — n 



i . . , . . i 



. , ' ^ ( 1 )' (35) 

In summary, for l-in-2 satisfiability on a chain, which is a simple XORSAT problem 19j with 
a tree factor graph, the belief propagation and susceptibility propagation give the correct 
magnetization and susceptibility. 

Consider now the same problem on the simplest graph with one loop, a ring. 

Namely, let G be a 1-dimensional ring, which is defined by identifying variable i — with 
N and adding a factor a = N— | as well as two incident edges (N, N—~) and (N — l, N—h). 
Moreover, we assume that iV is an even integer so that there is no frustration. 

BP has a continuous family of fixed points: 

n« !=fiW =(-l)M ± , (36) 

where A± is a constant [17J • As a consequence of the existence of this family of fixed points, 
(1 — M) is not invertible; in fact it has an eigenvector with zero eigenvalue, yo = (1,-1) 
where 1 corresponds to the //-block and — 1 corresponds to the fj block. In agreement with 
the existence of this dangerous eigenvector, one finds that the susceptibility propagation 
update rule does not converge. As the susceptibility messages are updated, r). i+ \ . picks 
up the constant shift 5ij = 1. This effect is accumulated as the messages go around the 
ring, and the consequence is that the messages diverge as t — > oo. 

In summary, for l-in-2 satisfiability on a ring, which is an XORSAT problem on a graph 
with a loop, the belief propagation can converge to a family of solutions for the magnetization 
among which only one solution is exact. On the other hand, the susceptibility propagation 
update does not have a fixed point, it diverges. In the simple case of a ring, this behaviour 
can be cured by using the finite temperature version of the BP and susceptibility propagation 
update equations. But in general there is no guarantee of convergence of loopy BP and loopy 
susceptibility propagation, and when they converge the quality of their results cannot be 
assessed a priori. FigJH gives an example of analysis of a small instance of l-in-4 satisfiability, 
giving an idea of the errors made by susceptibility propagation on small factor graphs. On 
the other hand, as for standard BP, one may hope that the method becomes better for large 
instances when the factor graph is locally tree-like. 
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FIG. 1: Comparison between the 2-point connected correlation function calculated exactly and that 
estimated with susceptibility propagation. A l-in-4 satisfiability instance on a randomly generated 
factor graph with N = 27 variables and M = 16 constraints with Poisson degree distribution with 
average degree I = 2.4856. 

IV. NUMERICAL INVESTIGATION OF SUSCEPTIBILITY PROPAGATION 

IN LOCKED OCCUPATION MODELS 

In this section we study the use of susceptibility propagation, together with decimation, 
in some locked occupation models. Specifically, we shall study random instances of a locked 
occupation problem, where the factor graph is uniformly chosen among the graphs with 
the following degree distribution. All function nodes have degree K and the variables have 
random degrees chosen from truncated Poisson degree distribution 



The basic message-passing algorithm that we use is described by the following pseudocode: 
Input : Factor graph, constraint- vector, convergence criterion, initial messages 




(37) 



for which the average degree is 




(38) 



Output : Estimate for 2-point connected correlation functions (or ERROR-NOT- 
CONVERGED) 



• Initialize messages 

• Repeat until everything converges 

— Update cavity fields and cavity biases vf^J^Xi) and ^i^(xi) with fl2D,d3D 

— Update cavity susceptibilities vf}+ a j{xi,Xj) and vfLi ji x i> x j) "with (JEJldHJ) with the 
help of vf\ a (xi) and i^Lii. x i) obtained above 

• Compute 1-variable marginals Pi(xi) from the fixed-point messages Pa->ii x i) by (BJ 

• Compute 2-point connected correlation functions p™ nn (xj, Xj) from the fixed-point mes- 
sages i>i%{xi) and v { *l^ a (xi,Xj) by (JT2J) 

This algorithm requires a memory proportional to kMN, and each step of iteration requires 
a computation of 0(N 2 ) for fixed k. 



A. Decay of correlations 

Figj2] shows the distribution of magnitude of 2-point connected correlation function com- 
puted with susceptibility propagation for all pairs of points in a graph for a fixed distance 
between the points. One observes a broad dispersion of correlations, and an approximate ex- 
ponential decay with the distance. Here we measure the distance d with the convention that 
each edge connecting a variable to a constraint is of length 1. Because of this exponential 
decay, it is possible to use in some cases approximate versions of susceptibility propagation 
which are faster and use less memory. This is done by truncating to zero the cavity suscep- 
tibilities Ui^aj.Oa^ij beyond some prescribed distance dist(a, j) > d or dist(z,j) > d and 
keeping only the correlation functions between pairs of variables not far from each other. 
Although one can estimate the 2-variable marginal distribution solely from the knowledge 
of cavity fields [IT, Sec. 14. 2. 3], this truncation provides us with a more efficient practical 
method to compute the 2-variable correlations between variables with d > 4. 
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FIG. 2: This graph shows how the 2-point connected correlation p|° nn (xj, Xj) decays as the distance 
d between xi and xj increases. At each distance d, the distribution of |p™ nn (0, 0)| is plotted. In 
the inset, logarithm of the average of that quantity is plotted against the distance. The instance 
is l-in-4 Satisfiability on a random factor graph with N = 1618 variables and M = 1000 factors 
with the truncated Poisson degree distribution with average degree I = 2.4856. 

B. Pair Decimation Algorithm 

As we mentioned in the introduction, decimation consists in finding a variable with the 
smallest entropy and fixing it to the most probable value. Assuming that the susceptibility 
propagation provides us with the good estimate for the 2-point connected correlation, we can 
think of decimation which acts on a pair of variable instead of a single variable. Let Xi and Xj 
be variables. If one defines a random variable y^ = l(xi = Xj), one can compute the entropy 
for yij once one knows the 2- variable marginal Pij(xi,Xj) = Pi{xi)pj(xj) + p^ nn (xi,Xj). In 
pair decimation, one identifies the pair with the smallest entropy of and one fixes 
either Xi = Xj or Xi + Xj = 1, depending on which event is the most probable according to 
the measured correlation. This results in a reduced smaller CSP, which is still an occupation 
problem. The efficiency of this novel decimation process depends on whether we can find a 
pair with less entropy than the single variable with the smallest entropy. It is easy to see 
that, in the absence of correlations, namely if Pij(xi, Xj) = pi(xi)pj(xj), then the entropy of 

is larger than the one of X{ or xj. So the whole procedure relies on being able to detect 
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FIG. 3: Comparison between the minimum entropy min(5'i, Sj) (where Si and Sj are the entropies 
of Xi, Xj )and S™, that of yij = \{xi = Xj). The instance is l-in-4 Satisfiability on a random 
factor graph with iV = 1618 variables and M = 1000 factors with the truncated Poisson degree 
distribution with average degree I = 2.4856. 

correlations. Figj3] shows that strongly correlated pairs can be found. 

In practice, we have used the following decimation algorithm which mixes the two strate- 
gies of single-variable decimation and pair decimation: 

Input : Factor graph, constraint- vector, convergence criterion, initial messages 

Output : A satisfying assignment (or FAIL-NOT-FOUND) 

• While graph has more than R variables: 

— Compute local entropy estimates for the 1-variable marginals 

— Compute local entropy estimates for the 2-variable marginals 

— if 'heuristic criterion finds that single-variable decimation is better', 

* then fix the value of the variable. 

* else identify a variable in the pair with the other (or its negation) 

— Locate completely locked nearest neighbor pairs 

— Clean the graph 



* Fix the value of isolated variables 

— Do warning propagation. 

— Identify local locked pairs 

• When the number of variables is equal to or smaller than R: perform an exhaustive 
search for satisfying assignments. If found 

— Then return the satisfying assignment 

— Else return FAIL-NOT-FOUND 

The heuristic criterion that we use in order to decide between the two types of decimation 
is the following. We locate a variable with the least entropy and a pair of variables with 
the least entropy for y^. When the former is less than 5th or is smaller than the latter, we 
choose to do single-variable decimation. 

For the optimal reduction of the entropy within a decimation step, it is reasonable to set 
5 t h = 0. However, we find that 5 t h > performs better for finding a satisfying assignment. 
The optimal value of 5 t h depends on the type of locked occupation model and the average 
degree. This fact can be interpreted as follows: the estimation of 1- variable marginals is more 
precise than the 2-variable ones within given computational resource, thus it is advantageous 
to respect the former if it is decisively small. 



Warning propagation is a message-passing algorithm described in |8|, |20|. It logically 
infers the value of variables one by one from local structure of the factor graph. 

In the identification of local locked pairs, we look at each degree-2 constraint and see if 
the constraint enforces = or y^ = 1. If it is the case, we identify this pair. 

The threshold for exhaustive search has been fixed in our simulations to R = 16. The 
performance of this algorithm is shown for l-in-4 satisfiability A = (0,1,0,0,0) and l-or-4-in- 
5 satisfiability A = (0, 1, 0, 0, 1, 0) in FigJH For l-in-4 satisfiability, data with randomization 
is presented: instead of fixing the most polarized variable or pair, we fix a variable or pair 
randomly chosen among a fixed number (here we adopt 8) of most polarized variables/pairs. 
The figure also shows the two important thresholds for these problems, which are values 
of the average degree (a measure of the number of constraints) separating qualitatively 
distinct phases. The probability that a satisfying assignment exists drops from 1 to 
at the 'satisfiability threshold' £ s in the large factor graph limit. Between the 'clustering 
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FIG. 4: Success probability of pair decimation process for l-in-4 satisfiability A = (0,1,0,0,0) 
(left) and l-or-4-in-5 satisfiability A = (0,1,0,0,1,0) (right) on a random factor graph with M 
constraints and average degree £ (right), plotted versus I . For comparison, the performance 
of simple belief-guided decimation process is shown. The vertical lines show the clustering and 
satisfiability thresholds. 

threshold' ( 6] and the satisfiability threshold, although the satisfying assignments still exist 
with probability one, it is very difficult to find one by the algorithms known so far, because of 
the splitting of the set of solutions into clusters. In both LOPs the performance is improved 
compared to the simple belief-guided decimation employed in 6|. Especially for l-or-4-in-5, 
the present algorithm works well above the clustering threshold, a region of I where all 
known algorithms are reported to perform poorly {(j). 



V. CONCLUSION AND DISCUSSION 



We have shown how to find satisfying assignments for locked occupation problems based 
on the measurement of correlation among variables. This is in contrast with the conventional 
method which is guided by 1-variable marginals only. Since flipping a variable in a LOP 
forces another variable far apart to be flipped, the performance of the algorithm is improved 
when we take the correlations into account. 

We have calculated correlations with the susceptibility propagation. In this method, the 
correlations between variables which ar efar apart can be calculated as well as between those 
which are neighbors. Namely, the convergence property is controlled by a single matrix Ai. 

The susceptibility propagation, however, requires more computational time and memory 



resource than the simple belief propagation. Therefore, as the problem becomes larger, we 
face a (polynomial) increase of computation time. The truncation introduced in subsec- 
tior JIV Al might give a remedy since it reduces by a factor of N the computation time as 
well as the memory use. The decay of correlation suggests that this is a reasonable ap- 
proximation. We have performed preliminary experiments to find the performance of this 
approximate algorithm. As expected, it behaves similarly to that without truncation, the 
performance being only slightly degraded. 
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